home *** CD-ROM | disk | FTP | other *** search
/ Atari Mega Archive 1 / Atari Mega Archive - Volume 1.iso / gnu / othergnu / bc102s.zoo / libmath.b < prev    next >
Text File  |  1992-03-04  |  5KB  |  256 lines
  1. /* libmath.b for bc for minix.  */
  2.  
  3. /*  This file is part of bc written for MINIX.
  4.     Copyright (C) 1991, 1992 Free Software Foundation, Inc.
  5.  
  6.     This program is free software; you can redistribute it and/or modify
  7.     it under the terms of the GNU General Public License as published by
  8.     the Free Software Foundation; either version 2 of the License , or
  9.     (at your option) any later version.
  10.  
  11.     This program is distributed in the hope that it will be useful,
  12.     but WITHOUT ANY WARRANTY; without even the implied warranty of
  13.     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14.     GNU General Public License for more details.
  15.  
  16.     You should have received a copy of the GNU General Public License
  17.     along with this program; see the file COPYING.  If not, write to
  18.     the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  19.  
  20.     You may contact the author by:
  21.        e-mail:  phil@cs.wwu.edu
  22.       us-mail:  Philip A. Nelson
  23.                 Computer Science Department, 9062
  24.                 Western Washington University
  25.                 Bellingham, WA 98226-9062
  26.        
  27. *************************************************************************/
  28.  
  29.  
  30. scale = 20
  31.  
  32. /* Uses the fact that e^x = (e^(x/2))^2
  33.    When x is small enough, we use the series:
  34.      e^x = 1 + x + x^2/2! + x^3/3! + ...
  35. */
  36.  
  37. define e(x) {
  38.   auto  a, d, e, f, i, m, v, z
  39.  
  40.   /* Check the sign of x. */
  41.   if (x<0) {
  42.     m = 1
  43.     x = -x
  44.   } 
  45.  
  46.   /* Precondition x. */
  47.   z = scale;
  48.   scale = 4 + z + .44*x;
  49.   while (x > 1) {
  50.     f += 1;
  51.     x /= 2;
  52.   }
  53.  
  54.   /* Initialize the variables. */
  55.   v = 1+x
  56.   a = x
  57.   d = 1
  58.  
  59.   for (i=2; 1; i++) {
  60.     e = (a *= x) / (d *= i)
  61.     if (e == 0) {
  62.       if (f>0) while (f--)  v = v*v;
  63.       scale = z
  64.       if (m) return (1/v);
  65.       return (v/1);
  66.     }
  67.     v += e
  68.   }
  69. }
  70.  
  71. /* Natural log. Uses the fact that ln(x^2) = 2*ln(x)
  72.     The series used is:
  73.        ln(x) = 2(a+a^3/3+a^5/5+...) where a=(x-1)/(x+1)
  74. */
  75.  
  76. define l(x) {
  77.   auto e, f, i, m, n, v, z
  78.  
  79.   /* return something for the special case. */
  80.   if (x <= 0) return (1 - 10^scale)
  81.  
  82.   /* Precondition x to make .5 < x < 2.0. */
  83.   z = scale;
  84.   scale += 4;
  85.   f = 2;
  86.   i=0
  87.   while (x >= 2) {  /* for large numbers */
  88.     f *= 2;
  89.     x = sqrt(x);
  90.   }
  91.   while (x <= .5) {  /* for small numbers */
  92.     f *= 2;
  93.     x = sqrt(x);
  94.   }
  95.  
  96.   /* Set up the loop. */
  97.   v = n = (x-1)/(x+1)
  98.   m = n*n
  99.  
  100.   /* Sum the series. */
  101.   for (i=3; 1; i+=2) {
  102.     e = (n *= m) / i
  103.     if (e == 0) {
  104.       v = f*v
  105.       scale = z
  106.       return (v/1)
  107.     }
  108.     v += e
  109.   }
  110. }
  111.  
  112. /* Sin(x)  uses the standard series:
  113.    sin(x) = x - x^3/3! + x^5/5! - x^7/7! ... */
  114.  
  115. define s(x) {
  116.   auto  e, i, m, n, s, v, z
  117.  
  118.   /* precondition x. */
  119.   z = scale 
  120.   scale = 1.1*z + 1;
  121.   v = a(1)
  122.   if (x < 0) {
  123.     m = 1;
  124.     x = -x;
  125.   }
  126.   scale = 0
  127.   n = (x / v + 2 )/4
  128.   x = x - 4*n*v
  129.   if (n%2) x = -x
  130.  
  131.   /* Do the loop. */
  132.   scale = z + 2;
  133.   v = e = x
  134.   s = -x*x
  135.   for (i=3; 1; i+=2) {
  136.     e *= s/(i*(i-1))
  137.     if (e == 0) {
  138.       scale = z
  139.       if (m) return (-v/1);
  140.       return (v/1);
  141.     }
  142.     v += e
  143.   }
  144. }
  145.  
  146. /* Cosine : cos(x) = sin(x+pi/2) */
  147. define c(x) {
  148.   auto v;
  149.   scale += 1;
  150.   v = s(x+a(1)*2);
  151.   scale -= 1;
  152.   return (v/1);
  153. }
  154.  
  155. /* Arctan: Using the formula:
  156.      atan(x) = atan(c) + atan((x-c)/(1+xc)) for a small c (.2 here)
  157.    For under .2, use the series:
  158.      atan(x) = x - x^3/3 + x^5/5 - x^7/7 + ...   */
  159.  
  160. define a(x) {
  161.   auto a, e, f, i, m, n, s, v, z
  162.  
  163.   /* Special case and for fast answers */
  164.   if (x==1) {
  165.     if (scale <= 25) return (.7853981633974483096156608/1)
  166.     if (scale <= 40) return (.7853981633974483096156608458198757210492/1)
  167.     if (scale <= 60) \
  168.       return (.785398163397448309615660845819875721049292349843776455243736/1)
  169.   }
  170.   if (x==.2) {
  171.     if (scale <= 25) return (.1973955598498807583700497/1)
  172.     if (scale <= 40) return (.1973955598498807583700497651947902934475/1)
  173.     if (scale <= 60) \
  174.       return (.197395559849880758370049765194790293447585103787852101517688/1)
  175.   }
  176.  
  177.   /* Negative x? */
  178.   if (x<0) {
  179.     m = 1;
  180.     x = -x;
  181.   }
  182.  
  183.   /* Save the scale. */
  184.   z = scale;
  185.  
  186.   /* Note: a and f are known to be zero due to being auto vars. */
  187.   /* Calculate atan of a known number. */ 
  188.   if (x > .2)  {
  189.     scale = z+4;
  190.     a = a(.2);
  191.   }
  192.    
  193.   /* Precondition x. */
  194.   scale = z+2;
  195.   while (x > .2) {
  196.     f += 1;
  197.     x = (x-.2) / (1+x*.2);
  198.   }
  199.  
  200.   /* Initialize the series. */
  201.   v = n = x;
  202.   s = -x*x;
  203.  
  204.   /* Calculate the series. */
  205.   for (i=3; 1; i+=2) {
  206.     e = (n *= s) / i;
  207.     if (e == 0) {
  208.       scale = z;
  209.       if (m) return ((f*a+v)/-1);
  210.       return ((f*a+v)/1);
  211.     }
  212.     v += e
  213.   }
  214. }
  215.  
  216.  
  217. /* Bessel function of integer order.  Uses the following:
  218.    j(-n,x) = (-1)^n*j(n,x) 
  219.    j(n,x) = x^n/(2^n*n!) * (1 - x^2/(2^2*1!*(n+1)) + x^4/(2^4*2!*(n+1)*(n+2))
  220.             - x^6/(2^6*3!*(n+1)*(n+2)*(n+3)) .... )
  221. */
  222. define j(n,x) {
  223.   auto a, d, e, f, i, m, s, v, z
  224.  
  225.   /* Make n an integer and check for negative n. */
  226.   z = scale;
  227.   scale = 0;
  228.   n = n/1;
  229.   if (n<0) {
  230.     n = -n;
  231.     if (n%2 == 1) m = 1;
  232.   }
  233.  
  234.   /* Compute the factor of x^n/(2^n*n!) */
  235.   f = 1;
  236.   for (i=2; i<=n; i++) f = f*i;
  237.   scale = 1.5*z;
  238.   f = x^n / 2^n / f;
  239.  
  240.   /* Initialize the loop .*/
  241.   v = e = 1;
  242.   s = -x*x/4
  243.   scale = 1.5*z
  244.  
  245.   /* The Loop.... */
  246.   for (i=1; 1; i++) {
  247.     e =  e * s / i / (n+i);
  248.     if (e == 0) {
  249.        scale = z
  250.        if (m) return (-f*v/1);
  251.        return (f*v/1);
  252.     }
  253.     v += e;
  254.   }
  255. }
  256.